program Create4Adb
    
    use netcdf
    use iasing_osse
    
    implicit none
    
    integer :: jo, jl
    
    integer :: dbid
    character(len=128) :: dbfile
    integer :: nObs, natm, nlev
    type(r4a_ids) :: db_ids
    type(r60_dat) :: db_dat
    integer, allocatable :: list(:)
    real(JPR4), dimension(43,7) :: rolvl
    
    integer :: ilon, ilat, idat
    real(JPR4) :: sectta, secttas
        
    include "checknc.interface"
    include "r60_initializer.interface"
    include "r60_reader.interface"
    include "interpolProf_f.interface"
    
    ! Reading selected
    open(1, file='selected', status='old', action='read')
    read(1,'(A)') dbfile
    call checknc(nf90_open(dbfile, NF90_NOWRITE, dbid))
    db_ids = r60_initializer(dbid)
    read(1,*) nObs
    if (nObs == 0) then 
        call checknc(nf90_inquire_dimension(dbid, &
            db_ids%atm_did, len=nObs))
        allocate(list(nObs))
        list = (/(jo,jo=1,nObs)/)
    else
        allocate(list(nObs))
        read(1,*) list
    endif
    close(1)
    
    
100 format(I8,3I10)
101 format(2F10.3,2G15.6)
102 format(12I6)
103 format(2F10.3,6G15.7)

    open(10, file="atm4aosse.dsf", action='write')
    ! Reading observations
    do jo = 1, nObs
        db_dat = r60_reader(dbid, db_ids, list(jo)) 
        rolvl(:,1) = interpolProf_f(NLM60, db_dat%dim%lev, &
            db_dat%vpr%tem, NLM4A, L434A)
        rolvl(:,2) = interpolProf_f(NLM60, db_dat%dim%lev, &
            db_dat%vpr%q, NLM4A, L434A)
        rolvl(:,3) = interpolProf_f(NLM60, db_dat%dim%lev, &
            db_dat%vpr%co2, NLM4A, L434A)
        rolvl(:,4) = interpolProf_f(NLM60, db_dat%dim%lev, &
            db_dat%vpr%o3, NLM4A, L434A)
        rolvl(:,5) = interpolProf_f(NLM60, db_dat%dim%lev, &
            db_dat%vpr%co, NLM4A, L434A)
        rolvl(:,6) = interpolProf_f(NLM60, db_dat%dim%lev, &
            db_dat%vpr%ch4, NLM4A, L434A)
        rolvl(:,7) = interpolProf_f(NLM60, db_dat%dim%lev, &
            db_dat%vpr%so2, NLM4A, L434A)
            
        idat = mod(db_dat%coo%dat,1000000)
        ilon = int(db_dat%coo%lon*100)
        ilat = int(db_dat%coo%lat*100)
        
        sectta  = 1/cos(db_dat%opa%iza*DEG2RAD)
        secttas = 1/cos(db_dat%opa%sza*DEG2RAD)
        
        write(10,100) list(jo), idat, ilon, ilat
        write(10,101) db_dat%opa%tsu, db_dat%opa%psu, sectta, secttas
        write(10,102) 40, 6, 1, 2, 3, 5, 6, 9
        do jl = 4, NLM4A
            write(10,103) L434A(jl), rolvl(jl,:)
        end do
        
        
    enddo
    deallocate(list)
    
end program
